{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# haddock3进行蛋白-聚糖对接流程演示\n",
    "\n",
    "相关链接:\n",
    "\n",
    "- [github代码/数据链接](https://github.com/haddocking/protein-glycans)\n",
    "- [list index out of range问题](https://github.com/haddocking/haddock3/issues/755)\n",
    "- [sample factor 400问题](https://github.com/haddocking/haddock3/commit/c54c4b6002b5a19b64f5af49bcf9e2c3b8825d31)\n",
    "\n",
    "代码库结构:\n",
    "\n",
    "- analysis:包含论文用到的所有分析和绘图脚本\n",
    "- data:论文研究产生的所有数据\n",
    "  - df_ss_all_data.tsv:每个结构和对接场景,计算top1,5,10,50,100,200结构中是否至少有一个结构达到高/中/可接受/勉强接受质量\n",
    "  - df_clt_all_data.tsv:每个结构和对接场景,计算top1,2,3,4,5,10,20,30,40,50个簇中最好的4个结构中是否至少有一个结构达到对应质量\n",
    "- pdb_files:每个复合物体系需要的所有pdb文件\n",
    "  - xxxx_r_b.pdb:蛋白结合态结构\n",
    "  - xxxx_l_b.pdb:聚糖结合态结构\n",
    "  - xxxx_r_u.pdb:蛋白游离态结构\n",
    "  - xxxx_l_u.pdb:聚糖游离态结构\n",
    "  - xxxx_l_ensemble.pdb:聚糖预采样的结构集合,由mdref模块生成\n",
    "  - xxxx_analysis.pdb:实验对照结构\n",
    "  - xxxx_ti-aa.tbl:约束文件,所有界面残基均为主动残基\n",
    "  - xxxx_tip-ap.tbl:约束文件,蛋白界面残基为主动残基,聚糖所有残基为被动残基\n",
    "- cfg_files:对接流程配置文件\n",
    "  - bound_default_ti-aa.toml:结合数据集对接流程,使用默认打分函数+ti-aa约束\n",
    "  - bound_vdw_ti-aa.toml:结合数据集对接流程,使用vdW打分函数+ti-aa约束\n",
    "  - bound_vdw_tip-ap.toml:结合数据集对接流程,使用vdW打分函数+tip-ap约束\n",
    "  - unbound_vdw_tip-ap_clust.toml:非结合数据集对接流程,使用vdW打分函数+tip-ap约束\n",
    "  - unbound_ens_vdw_tip-ap_clust.toml:集合结构对接,使用vdW打分函数+tip-ap约束\n",
    "  - mdref-glycans-sf400x16.toml:使用mdref进行聚糖预采样\n",
    "\n",
    "## 前置条件\n",
    "\n",
    "- 下载官方论文的代码库\n",
    "- 安装haddock3,如果以前安装过,注意更新(源代码关键参数已修改):`git pull origin main`.\n",
    "- 安装openbabel(用于处理pdb.gz文件)\n",
    "\n",
    "## haddock3:mdref对GLYCAM-Web生成的聚糖结构进行预采样\n",
    "\n",
    "以1C1L为例,`1C1L_l_u.pdb`为GLYCAM-Web生成的聚糖,`mdref-glycans-sf400x16.toml`为mdref的执行流程配置文件.\n",
    "\n",
    "具体流程为:\n",
    "\n",
    "1. topoaa:生成拓扑结构\n",
    "2. mdref:400次采样,生成400个结构\n",
    "3. rmsdmatrix:计算rmsd矩阵\n",
    "4. clustrmsd:按rmsd聚类,生成20个簇\n",
    "5. seletopclusts:每个簇选择1个结构(簇中心结构)\n",
    "\n",
    "预采样计算时长约35min(4核)/13min(25核),占用7.5GB内存(25核).\n",
    "\n",
    "关于mdref各参数的详细说明,参见haddock3:[mdref](https://www.bonvinlab.org/haddock3/modules/refinement/haddock.modules.refinement.mdref.html)文档.\n",
    "\n",
    "根据自己电脑的CPU核心数,可适当修改toml文件中的ncores,加速计算.默认使用4核心."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#通用设置\n",
      "\n",
      "mode = \"local\" #个人笔记本上使用local模式,源文件为batch模式\n",
      "queue = \"short\"\n",
      "queue_limit = 100\n",
      "concat = 5\n",
      "run_dir = \"run_mdref-glycans-sf400x16\"\n",
      "ncores = 25 # 使用的CPU核心数,根据自己的电脑修改\n",
      "\n",
      "# 提交给mdref的分子结构,路径从源文件略有调整\n",
      "molecules = [\n",
      "\t\"1C1L_l_u.pdb\"\n",
      "]\n",
      "\n",
      "[topoaa]\n",
      "\n",
      "[mdref]\n",
      "nfle1 = 1\n",
      "fle_sta_1_1 = 1\n",
      "fle_end_1_1 = 7\n",
      "sampling_factor = 400\n",
      "nemsteps = 200\n",
      "waterheatsteps = 100\n",
      "watersteps = 20000\n",
      "watercoolsteps = 8000\n",
      "\n",
      "[rmsdmatrix]\n",
      "resdic_B = [1,2,3,4,5,6,7] #论文研究的所有结构长度为2-7个单元,因此统一设置为1~7\n",
      "\n",
      "[clustrmsd]\n",
      "criterion = 'maxclust'\n",
      "linkage = 'average'\n",
      "n_clusters = 20\n",
      "\n",
      "[seletopclusts]\n",
      "top_models = 1\n",
      "top_cluster = 20\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%cat mdref-glycans-sf400x16.toml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 执行mdref流程\n",
    "!haddock3 mdref-glycans-sf400x16.toml\n",
    "# 日志太长,手动清除了笔记本内的日志"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0_topoaa  2_rmsdmatrix\t4_seletopclusts  data  traceback\n",
      "1_mdref   3_clustrmsd\tanalysis\t log\n"
     ]
    }
   ],
   "source": [
    "!ls run_mdref-glycans-sf400x16/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster_10_model_1.pdb.gz  cluster_18_model_1.pdb.gz  cluster_6_model_1.pdb.gz\n",
      "cluster_11_model_1.pdb.gz  cluster_19_model_1.pdb.gz  cluster_7_model_1.pdb.gz\n",
      "cluster_12_model_1.pdb.gz  cluster_1_model_1.pdb.gz   cluster_8_model_1.pdb.gz\n",
      "cluster_13_model_1.pdb.gz  cluster_20_model_1.pdb.gz  cluster_9_model_1.pdb.gz\n",
      "cluster_14_model_1.pdb.gz  cluster_2_model_1.pdb.gz   io.json\n",
      "cluster_15_model_1.pdb.gz  cluster_3_model_1.pdb.gz   params.cfg\n",
      "cluster_16_model_1.pdb.gz  cluster_4_model_1.pdb.gz   seletopclusts.txt\n",
      "cluster_17_model_1.pdb.gz  cluster_5_model_1.pdb.gz\n"
     ]
    }
   ],
   "source": [
    "!ls run_mdref-glycans-sf400x16/4_seletopclusts/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将生成run_mdref-glycans-sf400x16工作目录,每个模块的过程和结果文件均分别保存.其中`4_seletopclusts`中最终生成20个簇中的1个分子压缩文件.\n",
    "\n",
    "将其合并为结构集合文件(也可以从官方源代码数据中直接下载):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20 molecules converted\n"
     ]
    }
   ],
   "source": [
    "!obabel run_mdref-glycans-sf400x16/4_seletopclusts/*.pdb.gz -O 1C1L_l_ensemble.pdb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## haddock3:flexref对结构集合进行分子对接\n",
    "\n",
    "`1C1L_tip-ap.tbl`为约束文件,`unbound_ens_vdw_tip-ap_clust.toml`为分子对接流程配置文件.注意修改ncores参数.\n",
    "\n",
    "分子对接计算时长约23min(25核),内存占用约11GB(25核),磁盘占用约1GB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!haddock3 unbound_ens_vdw_tip-ap_clust.toml\n",
    "# 日志太长,手动清除了笔记本内的日志"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "00_topoaa      04_clustrmsd\t 08_caprieval\tanalysis\n",
      "01_rigidbody   05_seletopclusts  09_rmsdmatrix\tdata\n",
      "02_caprieval   06_caprieval\t 10_clustrmsd\tlog\n",
      "03_rmsdmatrix  07_flexref\t 11_caprieval\ttraceback\n"
     ]
    }
   ],
   "source": [
    "!ls run_unbound_ens_vdw_tip-ap_clust/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`analysis/11_caprieval`目录中`capri_ss.tsv`统计了top200结构的分数/rmsd等数据,`report.html`将对接相关的所有信息进行了绘图,`summary.tgz`中保存了对接结构的pdb文件,pdb中有对接结果的相关信息.\n",
    "\n",
    "如果`report.html`长时间不显示绘图,将`<script src=\"https://cdn.plot.ly/plotly-2.32.0.min.js\"></script>`替换为`<script src=\"haddock3所在环境的python包路径/lib/python3.9/site-packages/plotly/package_data/plotly.min.js\"></script>`(如`/home/regen/ana/envs/haddock3`).\n",
    "\n",
    "我的电脑上分数最佳的结构的IL-RMSD为0.849埃,所在簇平均IL-RMSD为0.582埃."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 清理生成的文件\n",
    "!rm -r run_*\n",
    "!rm 1C1L_l_ensemble.pdb"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cctbx",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
